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We explore the region of small sea quark masses below mps/mv ~ 0.5 in two-flavor QCD using a mean-field 
improved clover quark action and an RG-improved gauge action at a ~ 0.2 fm on 12^ X 24 and 16^ x 24 lattices. 
We find that instability of the standard BiCGStab algorithm at small quark masses can be mostly removed by 
the BiCGStab(DS-I/) algorithm, which employs L-th minimal residual polynomials with a dynamical selection of 
L. We also find singular spikes of AH in the HMC algorithm at moderate values of At. Nature of the spike is 
studied. We also study finite-size effects and chiral properties of meson masses. 



1. Introduction 

A rapid increase of the computational cost 
and instabilities in simulation algorithms have so 
far limited simulations of QCD with dynamical 
Wilson-type quarks to quark masses correspond- 
ing to mps/mv <;0.6 |l|], to be compared with 
the physical value niT^jmp — 0.18. This limita- 
tion of the range of quark masses causes sizable 
ambiguities and systematic errors in the extrap- 
olation to the physical point. A related problem 
is that the mass dependences predicted by chiral 
perturbation theory (ChPT) have not been con- 
firmed in full QCD simulations [|l[||J^. 

In this report, we explore the light quark mass 
region extending our previous systematic study 
of two-flavor QCD Through an improvement 
of simulation algorithms, sea quark masses corre- 
sponding to mps/mv = 0.6-0.4 are studied on 
coarse lattices with a ~ 0.2 fm. 

We adopt a renormalization-group improved 
gauge action and a meanfield-improved clover 
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Figure 1. True residual at {mps /mv)vai = 0.27 
as a function of L and Nuer measured on a con- 
figuration at {mps /mv)sea = 0.6. 



quark action At /? = 1.8 (a ~ 0.2fm), we 
employ 12^ x 24 and 16^ x 24 lattices with the 
spatial size 2.4 and 3.2 fm. Simulation parame- 
ters and present statistics are summarized in Ta- 
ble |l|. Measurements are done at every 5 trajec- 
tories generated by the HMC algorithm. 



2. Inversion of the quark matrix 

In Ref. [Q, the BiCGStab algorithm is adopted 
to invert the quark matrix. At small quark 
masses, the BiCGStab sometimes fails to con- 
verge. While the CG algorithm is guaranteed 
to converge, it is quite time-consuming. We find 
that an extension of BiCGStab to L-th order min- 
imal residual polynomials, the BiCGStab(i) al- 
gorithm is more stable. The conventional 
BiCGStab corresponds to the case L = 1. 

A larger L is expected to lead to a better 
convergence. In practice, however, too large 
L is time-consuming and also frequently intro- 
duces other instabilities. Figure |] shows the L- 
dependence of the residual on a test configura- 
tion. Optimum values of L depend on the simu- 
lation parameters. 

To avoid a tuning of L at each simulation point, 
we employ the BiCGStab(DS-L) algorithm |]. 
This is an improvement of BiCGStab(L) in which 
an optimum L is dynamically selected. We find 
that BiCGStab (DS-L) is much more robust than 
the original BiCGStab at small quark masses 
without much increase of the computer time. 

3. HMC updates 

Another difficulty at small quark masses is the 
appearance of instabilities in MD evolutions of 
the HMC algorithm. Figure || is an example 
of the time history of AH = Htriai — Hoid at 
mps/mv = 0.5. We observe that AH some- 
times shows huge values ("spikes"). A similar 
phenomenon has been reported in |Q. A conse- 
quence of the spikes is a distorted distribution of 
e-'^" at e-^" « (see the inset of Fig. |). While 
the trial configurations with large AH are au- 
tomatically rejected by the HMC algorithm, the 

Table 1 

Simulation parameters and number of trajectories 
at /3 = 1.8. 

{mps/mv)sea 0.6 0.5 0.4 

Ksea 0.14585 0.14660 0.14705 

Ntraj 12^ X 24 3500 3500 2100 

16^ X 24 1400 750 
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Figure 2. Time history of AH and the corre- 
sponding histogram of exp(— AiJ). 
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Figure 3. AH and the quark force at a spike. 

distortion of the e~^^ distribution can introduce 
additional systematic errors (within the statisti- 
cal errors) through the accept/reject process. 

We find that the frequency of spikes decreases 
as we decrease the MD step size dt. We have also 
checked that spikes do not violate the reversibil- 
ity and area-preservation required by the HMC 
algorithm. 

In order to clarify the origin of spikes, we plot 
AH{p{t),U{t+\dt)) and \\D-^{D'^)-^cj)\\ for the 
quark force as functions of the MD time i on a 
configuration with a spike, where p is the conju- 
gate momentum of U and cf) is the pseudo-fermion 
field. From Fig. ||, we find that AH jumps to a 
large value when the quark force becomes large. 
We confirm that the quark force remains small 
when dt is small enough to remove spikes. 

A mechanism to induce large AH from large 
force is described in H. In a simple harmonic 
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Figure 4. Pseudo-scalar and vector meson masses. 

oscillator system with the frequency a discrete 
approximation of the MD evolution becomes un- 
stable for dt > 2/llj and AH diverges exponen- 
tially with the MD time. This suggests that the 
discretized MD evolution of HMC with fixed dt 
may become unstable if an effective lu cx force 
in QCD exceeds a critical value. To study the ori- 
gin of large force, we are currently investigating 
the eigenvalues of the quark matrix. 

4. Meson masses 

Our preliminary results for the pseudo-scalar 
and vector meson masses are shown in Fig. ^, to- 
gether with the previous results at the same /? [Q . 
We note that the new data points are consistent 
with the previous data in the overlapping region, 
but show a deviation from the chiral extrapola- 
tion fit curves of the previous data, at small quark 
masses. The difference amounts to 3% (2.50-) 
in the chiral limit. Comparing the results from 
12^ X 24 and 16^ x 24 lattices, we find that the 
finite size effects are small at mps/mv ~ 0.6, 
while, at mps/mv — 0.5, my on f6'^ x 24 is 
slightly lower by 2% {1.5a). A higher statistics 
is needed to draw more definite conclusions. 

Using the data shown in Fig. |[ we test the 
PCAC relation, 

^^l + ^logy + ay,y^-^ (1) 
2Bomq Nf (47r/ps)^ 

where Bq and a are unknown parameters and 
C = 1 from ChPT. Fixing fps to U = 93 
MeV and adopting the AWl quark mass for m,. 



the data can be well fitted with x^/df — 0.05 
when we treat C as a free parameter. However, 
we obtain C = 0.06(2). Correspondingly, when 
we fix C = 1, the data cannot be described by 
Eq. (|l|) {x'^/df « 40). Therefore, our data down 
to mps/mv = 0.4 do not show the logarithmic 
curvature predicted by ChPT. 

5. Conclusions 

We explored the light quark region of QCD 
down to mps/mv = 0.4 using improved Wilson- 
type quarks. We found that the BiCGStab(DS-i) 
algorithm is robust at small quark masses. Spikes 
in AH may introduce an additional systematic 
error but can be suppressed by setting dt small. 
The origin of the huge quark force is still un- 
der investigation. Preliminary analyses on meson 
masses did not find agreement with the logarith- 
mic behavior expected from chiral perturbation 
theory to one loop order. 
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